%% PRACTICAL: AUTOMATED PREPROCESSING IN OPT
% 
% In this practical we will use the OHBA Preprocessing Tool (OPT) to 
% perform fully automated preprocessing. This will go through the following 
% steps:
%
%     1) Conversion of data into SPM format
%     2) Downsampling
%     3) High pass filtering
%     4) Automated AFRICA denoising
%     5) Epoching
%     6) Automated outlier trial rejection (using a Fieldtrip tool) 
%    
% We will work with a single subject's data from a button press experiment. 
% This can be downloaded from:
%
% <http://www.fmrib.ox.ac.uk/~woolrich/button_press_data.tar.gz>
%
% Note that this contains the fif file:
%
% fifs/loc_S02_sss1.fif
%
% Which has already been SSS Maxfiltered and downsampled to 250 Hz, and 
% which we will use as the input into this analysis.
%
% MWW & APB 2014


%% SETUP THE MATLAB PATHS
%
% Sets the Matlab paths to include OSL. Change these paths so that they 
% correspond to the setup on your computer. You will also need to ensure 
% that fieldtrip and spm are not in your matlab path (as they are included
% within OSL).
    
osldir = '/Users/abaker/Dropbox/Code/osl_latest/osl_svn/osl';    

addpath(osldir);
osl_startup(osldir);


%% SPECIFY DIRECTORIES FOR THIS ANALYSIS
%
% Change the datadir and workingdir variables to correspond to the correct
% directories.

% This is the directory containing the downloaded experimental data
datadir = '/Users/abaker/Scratch/button_press_data';
    
% This is the directory in which the analysis files will be stored
workingdir = datadir;

if ~isdir(workingdir)
    mkdir(workingdir);
end


%% SET UP THE LIST OF SUBJECTS AND THEIR STRUCTURAL SCANS FOR THE ANALYSIS
%
% Specify a list of the fif files, structual files (not applicable for this
% practical) and SPM files (which will be created). It is important to make
% sure that the order of these lists is consistent across sessions. 
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:

% fif_files{1}=[testdir '/fifs/loc_S01.fif']; 
% fif_files{2}=[testdir '/fifs/loc_S02.fif']; 
% etc...

% spm_files{1} = [workingdir '/spm_files/loc_S01.mat'];
% spm_files{2} = [workingdir '/spm_files/loc_S02.mat'];
% etc...

% structural_files{1} = [workingdir '/structurals/loc_S01.mat'];
% structural_files{2} = [workingdir '/structurals//loc_S02.mat'];
% etc...

%%%
clear raw_fif_files fif_files structural_files;

% List of the existing fif files for each subject
raw_fif_files{1} = [datadir '/fifs/loc_S02.fif']; 

% Our raw fif file has already been Maxfiltered:
fif_files{1} = [datadir '/fifs/loc_S02_sss1.fif'];

% List of the structural files for subject
structural_files{1}=[datadir '/structs/anat.nii']; % leave empty if no .nii structural file available


%% SETTING UP AN OPT ANALYSIS
%
% This sets up an OPT struct to pass to osl_check_opt, by setting 
% appropriate fields in the OPT struct. Note that some fields are mandatory
% while others are optional (and will be automatically set to their default
% values).
%
% The osl_check_opt.m function should be used to setup the settings for
% OPT. This function will check the settings, and will throw an error if 
% any required inputs are missing, and will fill other settings that are 
% not passed in with their default values. The OPT structure can then be 
% passed to osl_run_opt to do an OPT analysis.
%
% On the Matlab command line type "help osl_check_opt" to see what the 
% mandatory fields are.  Note that you MUST specify:
%
% opt.raw_fif_files: A list of the existing raw fif files for subjects 
%                    (need these if you want to do SSS Maxfiltering)
%
% OR:
%
% opt.input_files: A list of the base input (e.g. fif) files for input into 
%                  the SPM convert call
%
% AND:
%
% opt.datatype: Specifies the datatype; i.e. 'neuromag', 'ctf', 'eeg';
%
% For more information, see
% <https://sites.google.com/site/ohbaosl/preprocessing/opt-under-construction>

%%%
% *SPECIFY REQUIRED INPUTS*
%
% List of input files and data type

opt.epoch.time_range = [-1 2]; 
opt.epoch.trialdef(1).conditionlabel = 'StimLRespL';
opt.epoch.trialdef(1).eventtype      = 'STI101_down';
opt.epoch.trialdef(1).eventvalue     = 11;
opt.epoch.trialdef(2).conditionlabel = 'StimLRespR';
opt.epoch.trialdef(2).eventtype      = 'STI101_down';
opt.epoch.trialdef(2).eventvalue     = 16;
opt.epoch.trialdef(3).conditionlabel = 'StimRRespL';
opt.epoch.trialdef(3).eventtype      = 'STI101_down';
opt.epoch.trialdef(3).eventvalue     = 21;
opt.epoch.trialdef(4).conditionlabel = 'StimRRespR';
opt.epoch.trialdef(4).eventtype      = 'STI101_down';
opt.epoch.trialdef(4).eventvalue     = 26;
opt.epoch.trialdef(5).conditionlabel = 'RespLRespL';
opt.epoch.trialdef(5).eventtype      = 'STI101_down';
opt.epoch.trialdef(5).eventvalue     = 13;
opt.epoch.trialdef(6).conditionlabel = 'RespLRespR';
opt.epoch.trialdef(6).eventtype      = 'STI101_down';
opt.epoch.trialdef(6).eventvalue     = 19;
opt.epoch.trialdef(7).conditionlabel = 'RespRRespL';
opt.epoch.trialdef(7).eventtype      = 'STI101_down';
opt.epoch.trialdef(7).eventvalue     = 23;
opt.epoch.trialdef(8).conditionlabel = 'RespRRespR';
opt.epoch.trialdef(8).eventtype      = 'STI101_down';
opt.epoch.trialdef(8).eventvalue     = 29;

%%%
% Coregistration settings - we're not doing coregistration here, but 
% normally you would if you want to do subsequent analyses in source space 
opt.coreg.do = 0; 

%%%
% Outlier rejection settings
opt.outliers.do = 1;


%% CHECK OPT SETTINGS
%
% Now we can call osl_check_opt to check the fields we have specified and
% to fill in any remaining optional inputs

opt = osl_check_opt(opt);

% Display these settings
display(opt);

%%
% *LOOK AT OPT SUB-SETTINGS*
% The OPT structure contains a number of subfields containing the settings 
% for the relevant stages of the pipeline. Note that each of these has a 
% "do" flag (e.g. opt.downsample.do), which indicates whether that part of 
% the pipeline should be run or not.
%
% Take a look at opt.dirname. This will be the name of the directory (full 
% path) where the OPT will be stored, and is given a ".opt" extension. Note 
% that each OPT directory is associated with an OPT run - if you rerun with 
% the same opt.dirname then this will overwrite an old directory, and the 
% old OPT results will be lost. Hence, you should ensure that you change 
% opt.dirname for a new analysis, if you want to avoid overwriting an old 
% one! 
        
disp('Maxfilter settings:');
disp(opt.maxfilter);

disp('Downsampling settings:');
opt.downsample

disp('AFRICA settings:');
disp(opt.africa);

disp('Highpass filter settings:');
disp(opt.highpass);

disp('Epoching settings:');
disp(opt.epoch);

disp('Outlier settings:');
disp(opt.outliers);

disp('Coregistration settings:');
disp(opt.coreg);


%% RUNNING THE OPT ANALYSIS
%
% Run the OPT analysis using osl_run_opt
opt=osl_run_opt(opt);


%% VIEWING THE OPT RESULTS
%
% Running the OPT analysis will create an OPT output directory (whose name 
% is the name set in opt.dirname with a ".opt" suffix added). This contains
% all you need to access the results of the analysis. Note that you can 
% load these into Matlab using the call:

opt = osl_load_opt(opt.dirname);

%%% 
% In particular, the OPT object contains a sub-struct named results, 
% (i.e. opt.results), containing:
%
% .logfile             (a file containing the matlab output) 
%
% .report              (a file corresponding to a web page report with 
%                        diagnostic plots)
%
% .spm_files           (a list of SPM MEEG object files corresponding to 
%                      the continuous data (before epoching), e.g. to 
%                      pass into an OAT analysis)
%
% .spm_files_epoched   (a list of SPM MEEG object files corresponding to the epoched data, e.g. to pass into an OAT analysis)
%
%
% It is highly recommended that you always inspect both the 
% opt.results.logfile and opt.results.report, to ensure that OPT has run 
% successfully.


%% VIEWING THE OPT REPORT
%
% Open the web page report indicated in opt.results.report in a web 
% browser. This displays important diagnostic plots. At the top of the file 
% is a link to opt.results.logfile (a file containing the matlab output) - 
% check this for any errors or unusual warnings.
%
% Then there will be a list of session specific reports. Here we have only 
% preprocessed one session. Open this link in your web browser, or in the 
% Matlab web browser:

open(opt.results.report.html_fname)

%%%
% This brings up the diagnostic plots for session 1. There are a number of 
% things to look out for:
%
% *Maxfilter*: Normally, the first thing shown would be the results of running SSS 
% Maxfilter (and associated bad channel detection). Since we have not run 
% that here there are no diagnostic plots to show for this.
%
% *Histogram of events corrected for button presses:* Shows you the number 
% of triggers found for each event code - check that these correspond to 
% the expected number of triggers in your experimental setup.
%
%
% *AFRICA*
%
% *Mains artefacts:* IC sensor maps (for both sensor types), spectra, and 
% time courses detected as being due to 50 Hz mains noise by AFRICA - check 
% that these have sensible frequency spectra with a peak at 50 Hz
%
% *EOG and ECG artefacts:* IC sensor maps (for both sensor types), spectra, 
% and time courses detected as being due to EOG or ECG artefacts by AFRICA. 
% These have been found due to their IC time courses having high 
% correlation with the corresponding EOG and ECG channels in the data - 
% check that these have sensible time courses (at least for EOG) and 
% topographies (for both EOG and ECG) [you will learn this by experience].
% 
% *High Kurtosis artefacts:* IC sensor maps (for both sensor types), 
% spectra, and time courses detected as having very high kurtosis over time 
% by AFRICA. Very high kurtosis is caused by having very "peaked" 
% distributions, and are more likely to be due to non-neuronal artefacts - 
% check that these have appropriately "bizarre" time courses and 
% topographies.
%
% *Outlier Detection:* Histograms and scatterplots before and after outlier
% detection. The scatterplots show the channels/trial number versus the 
% metric (e.g. "std") as red crosses before rejection and green crosses 
% after rejection. Channels/trials to be retained are indicated by green
% circles.


